Monte Carlo studies of quantum and classical annealing on a double-well 



in 
o 
o 

(N 

o 

Q 

m 



CZ2 



c 

O 

o 



> 

O 

(N 

in 
o 

03 



o 
o 



X 



Lorenzo Stella, 1 Giuseppe E. Santoro, 1,2 and Erio Tosatti 1 ' 2 

1 International School for Advanced Studies (SISSA), 
and INFM Democritos National Simulation Center, Via Beirut 2-4, 1-34014 Trieste, Italy 
2 International Centre for Theoretical Physics (ICTP), P.O.Box 586, T34014 Trieste, Italy 

(Dated: February 2, 2008) 

We present results for a variety of Monte Carlo annealing approaches, both classical and quantum, 
benchmarked against one another for the textbook optimization exercise of a simple one-dimensional 
double-well. In classical (thermal) annealing, the dependence upon the move chosen in a Metropolis 
scheme is studied and correlated with the spectrum of the associated Markov transition matrix. In 
quantum annealing, the Path-Integral Monte Carlo approach is found to yield non-trivial sampling 
difficulties associated with the tunneling between the two wells. The choice of fictitious quantum 
kinetic energy is also addressed. We find that a "relativistic" kinetic energy form, leading to a higher 
probability of long real space jumps, can be considerably more effective than the standard one. 
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I. INTRODUCTION 

Solving difficult classical combinatorial problems by 
adiabatically controlling a fictitious quantum Hamilto- 
nian is a fascinating route to quantum computation 1 , 
which is alternative to the main-stream approach^. The 
quantum annealing (QA) strategy^* is precisely based 
on the adiabatic switching off, as a function of time, 
of strong fictitious quantum fluctuations, suitably intro- 
duced by an extra kinetic term in the otherwise classical 
Hamiltonian under consideration. 

A number of authors have by now applied QA to at- 
tack a variety of optimization problems. Some amount 
of success was obtained in several cases, notably in 
the folding of off-lattice polymer models 5,6 , in the 2D 
random fsing model annealingS^, in the Lennard-Jones 
clusters optimization^*^, and in the Traveling Salesman 
Problemii. 

Despite these successes - mostly obtained by Path- 
Integral Monte Carlo (PIMC) implementations of QA 
- there is no general theory addressing and even less 
predicting the performance of a QA algorithm. Differ- 
ent optimization problems are characterized by different 
"energy landscapes" of barriers and valleys about which 
little is known but which clearly influence the anneal- 
ing process^. Success of QA crucially depends on the 
efficiency with which the chosen kinetic energy and asso- 
ciated fictitious quantum fluctuations explores and influ- 
ences these effective energy landscape. This is a uncom- 
fortable situation, in view of the fact that it is a priori 
not obvious or guaranteed that a QA approach should 
do any better than, for instance, Classical Simulated An- 
nealing (CA). Indeed, for the interesting case of Boolean 
Satisfiability problems^ - more precisely, a prototypi- 
cal NP-complete problem such as 3- SAT - a recent study 
has shown that PIMC-QA performs definitely worse than 
simple CAil 

In a previous papes^, henceforth referred to as I, we 
compared classical and quantum annealing approaches 
in their performance to optimize the simplest potential 



landscape V(x), namely the one-dimensional potential 
double-well. Classical annealing was implemented by 
means of a Fokker-Planck (FP) equation with a time- 
dependent (decreasing) temperature T(t). Quantum an- 
nealing was performed through the Schrodinger's equa- 
tion, both in real and in imaginary time, describing a 
particle with a time-dependent (decreasing) inverse mass 
m~ 1 (t) in the double-well potential V(x). Given the 
textbook simplicity of the problem, it was possible to 
correlate in detail the features of the landscape (posi- 
tion and shape of valleys, barrier height) with the out- 
come of the annealing process. One saw for instance that 
Fokker-Plank CA is sensitive only to the ratio Ay/B of 
the energy splitting Ay between the two wells, or val- 
leys, and the barrier height B, whereas on the contrary 
Schrodinger QA is sensitive to the Landau-Zener tunnel- 
ing gap, and hence, among other things, to the barrier 
width. 

Unfortunately, such a direct approach is applicable, 
for practical purposes, only to insignificantly small sized 
optimization problems. To appreciate the difficulty, it is 
enough to consider that the number of possible configura- 
tions, and hence the Hilbert space size, of a small 32x32 
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ically large number which forbids any direct approach 
based on deterministic state evolution. Hard instances 
of typical optimization problems involve an even larger 
number of variables. Because of that, alternative strate- 
gies, based on Monte Carlo sampling, are mandatory. 
That notwithstanding, the various deterministic types of 
dynamics, Fokker-Planck, Schrodinger, or Monte Carlo 
dynamics, classical or quantum, are not at all equiva- 
lent. There is in particular no relationship between the 
(physical) time appearing in the Fokker-Planck or in the 
Schrodinger equation, and the corresponding classical or 
quantum Monte Carlo time-step. Furthermore, many dif- 
ferent types of Monte Carlo (MC) dynamics are possi- 
ble. While sampling the same equilibrium probability 
distribution and hence providing the same equilibrium 
averages^, they will generally perform very differently 
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in an out-of-equilibrium situation. In order to appreciate 
that, given the large variability of the possible MC out- 
comes, it is once again wise to concentrate attention on 
landscapes which are well under control. 

In this paper we implement and investigate differ- 
ent Monte Carlo annealing strategies, both classical and 
quantum, to the very same optimization problem - the 
one-dimensional double-well - which was studied in I by 
deterministic (Fokker-Planck and Schrodinger) anneal- 
ing approaches. We will first show (Sec. [H] that the 
choice of the proposed move in a Metropolis MC scheme 
strongly influences the annealing performance. Such a 
strong influence of the move is seen to be correlated with 
the instantaneous spectrum of the Markov transition ma- 
trix associated to the MC scheme, as shown in detail in 
Sec. lII Al Next, we will turn to Quantum MC, specifically 
to Path-Integral Monte Carlo (PIMC), and show how the 
required tunneling dynamics is highly non-trivial. The 
outcome of an allegedly state-of-the-art PIMC anneal- 
ing of a simple double-well can even be very unsatis- 
factory, due to sampling difficulties of instanton events 
(see Sees. IIIII IIII Al IIIIBfl . Finally, we will analyze the 
choice of kinetic energy in the QA Hamiltonian, by test- 
ing the effectiveness of an alternative relativistic kind of 
dispersion, which proves much more effective than the 
usual non-relativistic choice f Sec. IIII O . Section Hvl will 
finally contain a summary of the main results some dis- 
cussion. Technical details on the spectral analysis of the 
classical Markov process, and on the generalization of 
the bisection algorithm of interest for the PIMC study of 
Sec. IIII CI are relegated to the two Appendices. 



II. CLASSICAL MONTE CARLO ANNEALING 
OF A DOUBLE- WELL. 



Let us start with the double-well one-dimensional po- 
tential on which we will test the different classical and 
quantum Monte Carlo (MC) annealing schemes pre- 
sented in this paper. We assume V(x) to describe two 
generally asymmetric wells of the form: 



V(x) = 



Vq K o4 + +OX 

i 2 2 \2 

T t (x —a ) 



for x > 
for x < 



(1) 



where a+ > a_ (both positive), and the linear term ax, 
with a > 0, splits the degeneracy of the two minima. 
(The discontinuity in the second derivative at the origin 
is of no consequence in our discussion.) The optimal 
value of the potential is obviously E opt = V(x-). To 
linear order in the small parameter a, the two minima 
are located at x± = ±a± — aa±/(8Vb), and the splitting 
between the two minima is Ay = V[x+) — V(xJ) = 
a(a++a_), while the second derivative of the potential at 
the two minima, to lowest order in a, is given by V"(x = 
x±) = 8Vb/aj_. A "symmetric" double- well potential is 
recovered for a + = a_. 
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FIG. 1: Classical annealing. Comparison between Fokker- 
Planck (deterministic) annealing (triangles), and several dif- 
ferent Metropolis Monte Carlo annealings which differ only 
by the choice of the proposed move: Box, Gaussian, or 
Lorentzian (see text). The potential is a symmetric double- 
well of Eq. Qwith a_ = a+ = 1, Vo = 1, and a linear term 
a — 0.1. The initial condition was To = Vo = 1. The solid 
line through the points is a fit with a power-law of leading ex- 
ponent T ~ A v/ B times logarithmic corrections, see I. Clearly, 
the Fokker-Planck and the Monte Carlo annealing "times" are 
not related (see text). 



Starting with some initial state or distribution, the key 
quantity under inspection will be the average residual en- 
ergy e res (t) — Epot (t) — E op t after annealing for a certain 
time r, where E pot (T) is the average final potential en- 
ergy. In the classical MC case, annealing is performed 
by reducing the temperature T(t) = T (l — t/r). As 
shown in I fRef . Ilfil) . the results of a Fokker-Planck an- 
nealing are roughly independent of the asymmetry of the 
potential, and only controlled by the ratio Ay jB of the 
splitting Ay between the two minima and the barrier 
height B = Vo — V(x + ). In Fig. we present the results 
of a standard Metropolis MC classical annealings for a 
symmetric (a_ — a + — 1, Vq = 1) double- well potential, 
where a linear term ax with a = 0.1 provides a splitting 
between the two minima of order Ay ~ 0.2. Entirely 
similar results (not shown) are obtained for asymmetric 
double- wells. The solid triangles represent the results of 
exact integration of Fokker-Planck's (FP) equation, via a 
fourth-order Runge-Kutta method, with a linear anneal- 
ing schedule for the temperature T(t) = T (l — t/r), as 
reported in I. The circles, squares, and diamonds repre- 
sent the results of three different classical Metropolis MC 
annealings - clearly showing very different behaviors -, 
differing uniquely in the choice of the proposed move. 
To clarify this point, we recall that the Markov process 
associated with a Metropolis MC is fully specified by a 
transition probability VF(a;'|a;), to move from x to x' in 
the time step At. W is the product of two factors: 



W{x'\x) = A(x'\x)P(x'\x) fon'^i 



(2) 
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where P(x'\x) is the probability of proposing a move from 
i to s', while 

A(x'\x)^Mm\l,^f^e-^- E ^ k - T ) , (3) 
I P\X \x) J 

is the probability of accepting the proposed move. The 
various results shown in Fig. ^ differ in the choice of 
P(x'\x), i.e.: 

1) (Box, circles) the new proposed position x' is uni- 
formly distributed inside a box of length 2a centered 
around x, i.e., 

P(x'\x) = P Bo x(x'\x) = -U(<7 - \x' - x\) ; (4) 
la 

2) (Gaussian, squares) the new position x 1 is distributed 
according to a Gaussian of width a centered around x, 
i.e., 



P(x'\x) =P Gau {x'\x) = 



1 



(5) 



3) (Lorentzian, diamonds) the new position x' is dis- 
tributed according to a Lorentzian of width a centered 
around x, i.e., 



Evidently, both the Gaussian and the Lorentzian trial 
moves allow the attainment of residual energies e res (r) 
much below that of the Box move (and of the FP dy- 
namics). The Lorentzian case is particularly impressive. 
If we fit the asymptotic part of the residual energy data 
by c res (r) oc r~ n ° A , the annealing exponent Vic a turns 
out to be ~ 0.31 for the Gaussian case, and ~ 1.0 for the 
Lorentzian one, i.e. definitely much larger than what the 
Huse and Fisher theory^ predicts for the Fokker-Planck 
case, namely Vic a — Ay/B ~ 0.2. 

The long Lorentzian tails clearly imply an abundance 
of moves operating long jumps, which are benefi cial to 
the annealing even at small temperatures T. Ref. Il8ll9l 
indeed have reported a genuine improvement of the sam- 
pling by making use of the Lorentzian distribution, in- 
stead of the Gaussian, for the Thompson problem (find- 
ing the equilibrium distribution of N charged particles on 
a sphere and in optimizing the structure of a cluster 
made of N Ni atomsi^. In order to analyze these issues in 
more quantitative detail, we carried out a spectral anal- 
ysis of the transition matrices W(a;'|x) corresponding to 
the different proposal move schemes. In the following 
section we briefly recall the crucial steps behind this ap- 
proach. 



P(x'\x) =Pl oi {x'\x) = 



1 



7T (x' 



(6) 



In all three cases, the value of a, which controls the 
move range, was decreased with temperature according 
to a = <7o \/T/To, with ao — 2.0, for the Box and the 
Gaussian case, and a$ = 2.9 for the Lorentzian one. 
This choice was made empirically by monitoring the av- 
erage acceptance and the root mean square displacement, 
y/ ({x i+ i — Xi) 2 ), of the MC process. 

As evident from Fig. the results of the various 
choices of proposal moves differ in a substantial way. 
While the Box choice leads to a residual energy very close 
to the deterministic Fokker-Planck result, the Gaussian 
choice, and even more so the Lorentzian one, lead to 
residual energies which decrease faster with the anneal- 
ing time r. 

The reason for the non uniqueness of the result of a 
classical MC annealing is not difficult to grasp. The 
Metropolis algorithm guarantees that the expectation 
values of a system at equilibrium are correctly reproduced 
by a MC time average, after a suitable equilibration time 
T eq . The equilibration time T eq depends on the quantity 
one is interested, and, more crucially for our scope, on the 
type of proposed moves, i.e., on the P{x'\x). For equi- 
librium MC simulations, this issue is of little concern: 
As long as one checks that the stochastic process has 
reached a stationary condition, the results are guaran- 
teed to be correct. In an annealing simulations, however, 
(or, in that respect, in any out-of-equilibrium evolution), 
the choice of the P{x'\x), determining the instantaneous 
equilibration time of the process, strongly influences the 
overall performance of the algorithm. 



Spectral analysis of the Markov process in 
Monte Carlo 



Every discrete-time Markov process can be defined by 
means of an equation of the form: 



P(i,t + At) = J2 W(i\j)P(j,t) 



(7) 



where, for simplicity of notation, we have assumed that 
the possible states of the system, labeled i and j, are also 
discrete (we will indeed reduce our continuous problem to 
a discrete one in the following numerical calculations, see 
below). P(i,t) is here the probability of finding the sys- 
tem in state i at time t, and W(i|j) is the so-called Tran- 
sition Matrix, i.e., the transition probability for going 
from state j to state i during the time step At. The tran- 
sition matrix M^(«|j) (assumed to be time-independent) 
is a stochastic matrb*2!l, in the sense that j) > and 
Si = 1. In our particular case, the stochastic pro- 

cess behind the Metropolis sampling is given by Eq. [3 so 
that W^(«|j) also fulfills the detailed balance condition 16 . 
The idea underlying the spectral analysis of is 
very simple. For a fixed temperature T, the Markov pro- 
cess admits a spectral analysis in much the same spirit in 
which one analyzes a time-dependent Schrodinger equa- 
tion via the spectrum of the associated Hamiltonian. In- 
deed, in continuous time - i.e., for At — ► - the result 
is quite simple to state. A similar analysis can be car- 
ried out for the original discrete-time Markov chain, see 
Ref. I20L As explained there, from W(i|j) one can define 
a closely related matrix W^(i|j) appearing in the short- 



time expansion, W^i|j) 



AtW(i\j), of the expo- 
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nential map: W = e Atw ; in turn, — W possesses a set 
of right eigenvectors {P k {i)} with P (i) = P^(i) (the 
equilibrium distribution), and the corresponding eigen- 
values {A*;} are such that A = 0, X k > Vfc > 0. In 
terms of the spectrum of W, one can formally write the 
general solution of the time-continuum (At — > 0) limit of 
Eq. © as^: 



P(i, t) = P^ (i) + e' lk 4 P k (*) 



(8) 



fe>0 



where a k are real coefficients depending on the initial 
(t = 0) condition. Eq. shows that the approach to the 
equilibrium distribution is governed by a relaxation time 
t re i oc A J" , Ai being the smallest non-zero eigenvalue of 
—W (recall that Ao = 0)_. Therefore, the larger the spec- 
tral gap A = Ai — Ao = Ai, the faster the system reaches 
equilibrium. However, since in the Metropolis algorithm 
the time step is actually discrete (see appendix IX}) . it is 
better to consider the set of the W eigenvalues, given 
by Ai = e _AtAi , with Ao = 1. Therefore, in the rest 
of the paper we shall refer to the spectral gap of the 
original transition matrix W, which is A — 1 — Ai, in- 
stead of the corresponding time-continuum limit intro- 
duced above. As we shall show in a while, this choice 
will not change the conclusion of our spectral analysis. 

How can we use the spectral gap concept in an anneal- 
ing context? It is natural to think that, for a given sched- 
ule T(t), it would be good to maximize the spectral gap 
for each given temperature T(t). This is indeed what we 
will explore in the following. Generally speaking, find- 
ing explicitly the spectral gap A for a given W is an 
impossible task; in the present one-dimensional double- 
well context, however, it turns out that the spectrum 
of W(a;'|a;) is quite simple to obtain by, for instance, a 
straightforward discretization of the real axis. 

We present now the results obtained by applying this 
spectral analysis to the proposal moves we defined in 
Sec. [H] We recall that they are: The Box move, see 
Eq. 01 the Gaussian move, see Eq. \5\ and the Lorentzian 
one, see Eq.|SJ The main quantity under inspection is the 
spectral gap A of the transition matrix W(x'|x) which is 
a function of both the temperature T and the proposal 
move range a, A(T, a). We obtain this gap (as well as 
the higher eigenvalues) by diagonalizing an appropriate 
symmetrized (see appendix 0) and discretized version 
of the transition operator W(x'\x)£* In Fig. |2 we plot 
the spectral gap A versus a at two fixed temperatures, 
T/Vq = 0.1 (Top) and T/V = 0.001 (Bottom). As a 
general remark, we note how the behavior of A(er) is re- 
markably different for different trial moves. Interestingly, 
for every proposal move and every T, there is a value of 
a which maximizes the gap A. Recalling that the re- 
laxation time is proportional to A" 1 we anticipate that, 
at every temperature T, there is a <j(T) which provides 
the fastest relaxation. Clearly this a(T) is a priori a 
better choice than taking just a(T) = oo (T/Tq) 1 / 2 , as 
previously implemented in Sec.[n] At low T however the 
dependence of A on a is no longer smooth. In particu- 
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FIG. 2: Exact diagonalization of the transition operator W , 
for the symmetric potential case. We plot the spectral gap 
A versus the proposal move range a. The temperature is 
T/Vo = 0.1 (Top) or T/Vo = 0.001 (Bottom), in the same 
units in which the barrier parameter is Vq = 1. The solid 
and dashed black lines represent the fit provided by Eq. IA8I 
(Gaussian case) and Eq . I A9I (Lorentzian case). The functional 
form of those fits is recalled in the legend. 



lar, for T/Vo — 0.001, the Box trial move shows a sharp 
transition from vanishingly small gap values to finite 
ones, while both Gaussian and Lorentzian moves show 
a cusp maximum. At low temperatures, the maximum 
gap AmaP of the Lorentzian move is larger than Amax^ , 
which is in turn larger than A^f Q °f'; correspondingly, for 
small a, we have A( Lor V) > M Gau \a) > A^ Box \a). 
This is, essentially, the reason for the faster relaxation 
dynamics of the Lorentzian move during annealing. No- 
tice also that the maximum value A max of A, decreases 
with temperature in all cases, as one can see by com- 
paring the two panels in Fig. |2 Moreover the value of 
a that maximizes the Lorentzian and Gaussian gap de- 
creases with temperature, while it seems to converge to 
a finite value (around a w 1.75) for the Box case. 

The Box case is particularly intriguing, and deserves a 
closer inspection. In Fig.|21we plot A versus inverse tem- 
perature 1 /T for the case of the Box move, and for several 
values of a. For a whole range of a, a < a cr ~ 1.75, we 
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FIG. 3: Exact diagonalization of the transition operator W , 
for the symmetric potential, with Box proposal move. We 
plot the spectral gap A versus the inverse temperature 1/T 
for several fixed values of the proposal move range a. Top 
left inset: The value of the effective barrier, B c s, seen by the 
system when a Box proposal move with range a is employed. 
The value a cr - which corresponds to a vanishing effective 
barrier - is indicated by a vertical arrow. Top right inset: 
The two first gaps, Ai = 1 — Ai and A2 = 1 — A2, for a Box 
proposal move versus the inverse temperature, 1/T, at fixed 
move range a = 1.78 



see a low-temperature behavior typical of an activated 
Arrhenius process A(T) oc e^ Bet( ^"'' T , where B c s(a) is 
an effective barrier, clearly a-dependent, which the sys- 
tem experiences. Above a certain a cr , the behavior of 
A changes drastically, from Arrhenius to what appears, 
at first sight, just a constant. A closer inspection on an 
extended temperature range, see Fig. (Top right in- 
set), shows that at very small temperatures an (avoided) 
crossing between eigenvalues of W has taken place, and 
A starts to decrease again toward zero as A oc T 1 / 2 ^ 
The effective barrier B c r((j) extracted from an Arrhenius 
fit, A oc e- B <>tt{a)/ T ; g 0es to zero for a — > a cr - There is 
a simple geometrical explanation for this critical behav- 
ior. First of all, we note that the value of a cr is close 
to the distance between the two minima. This suggests 
that the transition is related to the availability of a di- 
rect jump from the bottom of the metastable well to the 
other one. If we call x±, respectively, the minimum in 
the right and left well, there is a non-trivial solution of 
the equation V(x\) — V(x + ) with X\ lying between the 
two minima, x_ < x\ < x + (see Fig. 0J: It turns out 
that a cr = x+ — x%. 
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FIG. 4: Illustrative sketch of a Box move in a double-well po- 
tential. In the upper panel we consider the case a — a cr , while 
in the lower one we illustrate the general case a < a cr . Here 
V+ — V(x+) is the potential at the bottom of the metastable 
minimum. The meaning of the other symbols is explained in 
the text. 

Indeed, for a > a cr there is a possible proposed move 
that brings the system from x = x + (the metastable min- 
imum) to the point x' — x\ on the other side of the bar- 
rier: Such a "non-local" move pays no energy (AE = 0), 
and is therefore certainly accepted by the Metropolis al- 
gorithm. In some sense, for a > a cr the barrier is not 
seen (or, at least, there are allowed and accepted moves 
that do not see it), and B c s is zero. Consider now the 
case a < a cr . For every value V = V(x+) + AV, with 
AV > 0, there are two equipotential solutions xi(AV) 
and X2(AV), such that V(xi) — V{x-2) — V, lying be- 
tween the two minima, but on opposite sides of the bar- 
rier, i.e., with X- < xi(AV) < X2{AV) < x+. Denote 
now by d(AV) = x 2 (AV) - Xi(AV) the distance be- 
tween such two equipotential points (d{AV) is a mono- 
tonically decreasing function of AV). If the box width 
is cr, we can always find AV such that a = d(AV), i.e., 
such that there is a proposed move connecting x = x 2 to 
x' = x\ which, once again, bypasses the top of the bar- 
rier. The effective barrier seen by the system close to the 
metastable minimum, for a given value of cr, is therefore 
just (see Fig. @J 

B cS {a) = AV = d- 1 (a) , (9) 

i.e., in essence, a piece of length a is cut from the top of 
the barrier, and is effectively not seen by the system. In 
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other words, the effective barrier B c ff is just the potential 
energy drop AV^ which the system must overcome before 
a long jump | re' — ac ] ~ a is made available at no en- 
ergy cost (V(x') ~ V(x)). The value of B e g(a) obtained 
through such a simple geometric construction is shown 
by a solid line in Fig. Top left inset): The agreement 
with the numerical data - extracted from the Arrhenius 
fit - is remarkably good, with small deviations for very 
small values of a, which are likely due to the effect of 
the finite grid employed in diagonalizing the transition 
operator W (see Sec. Ill A|) . and to finite T effects»2i 

We would like stress that this geometrical picture is 
strictly true only at zero temperature, and does not apply 
to the Gaussian or Lorentzian cases, whose tails provide 
a small non- vanishing chance of making a "long jump" 
for any value of a. The treatment of the Gaussian and 
Lorentzian cases needs, therefore, a more specific and 
technical discussion, which we sketch in appendix lAl 

Summarizing, the Box trial move shows a sharp (first- 
order- like) transition at a value of o~ cr ~ 1.75, where the 
effective Arrhenius barrier B c ff(a) vanishes, and the gap 
starts to decrease as a power- law, A cx T 1 / 2 , for very 
small values of T (see Fig. |3J). 

One might suspect that the cusps shown in Fig. [21 for 
the Gaussian and Lorentzian cases signal, similarly to 
the Box case, some kind of transition. A closer inspec- 
tion, however, shows that this is not so: It is just a 
level crossing phenomenon. In Fig. [5] we plot A versus 
the inverse temperature 1/T for several values of a, for 
the Gaussian (Top panel) and Lorentzian (Bottom panel) 
proposal moves. After an initial Arrhenius-like behavior, 
particularly visible for the Gaussian small a cases, the 
system always, and smoothly, changes to a low-T behav- 
ior which is entirely similar to the large-a Box case: An 
apparent saturation of A followed by an avoided cross- 
ing between the first two exited states (not shown) and 
a final A cx T 1 / 2 . 22 No real transition exist in the Gaus- 
sian and Lorentzian case: The cusp in Fig. Amoves down 
toward smaller and smaller values of a for decreasing T. 
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FIG. 5: Exact diagonalization of the transition operator, for 
the symmetric potential, with Gaussian (Top) or Lorentzian 
(Bottom) move. We plot the gap A versus the inverse tem- 
perature 1/T for several fixed values of the proposal move 
range a. 



B. Classical Annealing with optimal a(T) 

In the previous section we argued about the possibility 
of an optimal choice for the proposal move range cr(T). 
This choice should guarantee the fastest instantaneous 
relaxation toward the instantaneous equilibrium distribu- 
tion P eq (x) at any given temperature. The classical an- 
nealing performance is expected to be greatly improved 
by such a choice of c(T). 

In practice, for each choice of proposal move (Box, 
Gaussian, Lorentzian), we numerically determined the 
value of a that maximizes the spectral gap A for each 
fixed temperature T. We will refer here to such an op- 
timal a as a op t (T) . In Fig. [5] we plot such an optimal 
choice for all the types of proposal moves introduced in 
Sec. [n] For comparison, we also show the o~(T) we em- 
ployed in Sec. [HJ which is definitely smaller than o~ op t, 



for given value of T. 

Having done this, we have then performed classical 
MC annealing runs, similarly to those reported in Sec.lTTl 
with the usual linear annealing schedule for the temper- 
ature T(t) = Tq (1 — t/r), where t is the annealing time. 
In Fig. 0we show the MC annealing results. We stress 
again that the difference with respect to the runs illus- 
trated in Fig.^ where we took a = ao(T ' /Tq) 1 / 2 , is that, 
here, o~ op t(T) is the optimal choice of a obtained from 
the maximum instantaneous spectral gap of the transi- 
tion matrix. We notice that the Box and Gaussian re- 
sults are now very different form the previous ones (see 
Fig. Moreover, interestingly, Box and Gaussian data 
fall almost on the same curve, which is in turn very close 
to the Lorentzian case, the latter showing once again the 
best annealing results^ 
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FIG. 6: Plot of the optimal a op t(T) for the Box, Gaussian 
and Lorentzian proposal move. We also show, for comparison, 
the schedule a(T) cx T 1/2 employed in Sec. |TT| 
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FIG. 7: Plot of Monte Carlo classical annealings for the sym- 
metric potential, r is the annealing time and e re3 the residual 
energy (see text). We report the results of the exact integra- 
tion (Fokker-Planck), together with the actual MC data for 
several proposal moves (Box, Gaussian, Lorentzian) . The MC 
simulations are performed with an optimal choice for er(T) ob- 
tained from the maximum instantaneous gap (see text). 



III. QUANTUM ANNEALING: 
PATH-INTEGRAL MONTE CARLO OF THE 
DOUBLE- WELL 



In I (Ref. 1 151) we studied, by direct numerical integra- 
tion, the real- and imaginary-time Schrodinger dynamics 
provided by the time- dependent Hamiltonian 



H 



-r(t) 



dx 2 



V(x), 



(10) 



where T(t) = ti 2 / [2m(t)] = To (l — -) is the inverse mass 
parameter appearing in the usual non-relativistic kinetic 
energy, providing the strength of the quantum fluctua- 
tions. This is a literal implementation of the Quantum 



Annealing (QA) strategy. However as pointed out in the 
introduction, this Schrodinger dynamics is suitable only 
for toy problems with a very limited Hilbert space, and 
quite inapplicable for actual optimization. In order to be 
a viable strategy for realistic optimization, QA must re- 
sort to stochastic, i.e., Quantum Monte Carlo (QMC), 
approaches, mostly appropriate for an imaginary-time 
framework (we remind the reader that, as shown in I, 
working in imaginary-time is actually beneficial for Q A) . 

There are several possible QMC techniques on which 
a QA strategy can be build. By far, the simplest of 
these QMC techniques is the Path-Integral Monte Carlo 
(PIMC) method, which has already been used with 
some success m the QA context 7 ^ 11 ! 14 . The method 
does not addresses the imaginary- time Schrodinger time- 
dependent dynamics, but simulates an equilibrium quan- 
tum system, held at a small finite temperature T, 
where the relevant quantum parameter T(t) is externally 
switched off in the course of a QMC simulation: t is not 
treated, therefore, as a proper physical time, and is only 
replaced by a Monte Carlo timem^ 

In the present section we will explore the potential 
of a PIMC-based QA strategy on the same simple one- 
dimensional double- well potential. The reason for invest- 
ing so much effort on this simple problem, is that the 
simplicity of the problem landscape will allow us to have 
a perfect control of all the ingredients of the method. In- 
deed, we will learn a lot about PIMC-QA, particularly 
about the limitations of the method, from investigating 
this toy problem. 

The equilibrium properties of our quantum system, for 
any fixed value of F = T(t), are all encoded in the quan- 
tum partition function 



Zifi) = Tr 



-/3(T+V) 



where 



T= -T 



dx (x\e-^ f+ir) \x) (11) 



(12) 



dx 2 



from which all the thermodynamics follows. Here, as 
usual, f3 = (Ub T) -1 . We will drop the Boltzmann con- 
stant ks from now on. As we see from Ea. llll Z involves 
an integral over a positive distribution, {x\e~^ T+v > \x), 
which involves, however, the very difficult task of calcu- 
lating the diagonal matrix element of the exponential of 
the Hamiltonian H. 

The standard approach is to rewrite Ea. llll as a Feyn- 
man Path-Integral . The main mathematical tool em- 
ployed in such an approach is the so-called Trotter theo- 
rem, which reads: 



(13) 



where P is the number of Trotter's slices (or replicas), in 
which the imaginary time interval [0, h/3] has been par- 
titioned, each slice being of length At = (H/3)/P. By 
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means of Eq. and inserting P — 1 identities in the 
form of li = J dxi\xi)(xi\, one can rewrite the partition 
function (see, for instance, Ref. \2l\ for this simple deriva- 
tion) as: 

Z(ff) = (j^) 2 / if d Xi e- s ^+0(At 2 (3) , 

(14) 

where Spa[x] is the so-called (euclidean) primitive ac- 
tion: 




with periodic boundary conditions xp = Xo, as a conse- 
quence of the trace present in Ea.HTl) 

The similarity of Eq. [21 with the classical partition 
function of a closed polymer with P beads is evident: 
The polymer beads Xi can be seen as the imaginary-time 
snapshots Xi = x(iAt), i = 0, • • • , P — 1, of a fluctuat- 
ing closed path x(t) in the enlarged configuration space 
(x,t), where t is the imaginary time. Two neighboring 
beads x\ and X{-\-% interact with harmonic interactions of 
spring constant K 1 - = mP 2 /(h(3) 2 , originating from the 
propagator of the quantum kinetic term, 

(xi\e-$ f \x l+1 ) oc e-^^+i"^) 2 . (16) 

The strength of the harmonic interactions in imaginary- 
time controls the amount of quantum fluctuations: A 
large mass m (classical regime) results in a strong K 1 - 
and hence in a very "rigid" polymer, while a small mass 
m (quantum regime) results in a very soft and fluctuating 
polymer. All beads are also subject to the classical po- 
tential V{x). Once we have reduced our problem to the 
Path-Integral form in Eq. [21 the standard techniques of 
classical Metropolis MC can be used, and the resulting 
algorithm is what is called a Path-Integral Monte Carlo 
(PIMC). The most obvious MC moves to be used are just 
single bead moves, exactly as in a classical MC. 

These are the bare bones of a PIMC approach. For 
the problem we are dealing with, a particle in a poten- 
tial (or, more generally, for systems of quantum particles 
on the continuum), one can improve on the method just 
sketched in two possible directions: i) by improving the 
quality of the approximation in Eq. 1141 so as to get a 
smaller Trotter discretization error— for a given number 
of Trotter slices P used^S U) by introducing MC moves 
which arc more sophisticated than just moving a single 
bead at a time. Regarding i), we will adopt a fourth or- 
der approximation to the action^S^i, which improves the 
Trotter truncation error of the partition function from 
O (A t 2 0) to O (A t 4 0) . More precisely, we used the 
so-called Takahashi-Imada approximation 30 ' 3 ^, which is 
especially suitable for continuous systems. The partition 
function is still given by an expression of the form Eq. 1141 
with the primitive action Spa replaced by the following 



Takahashi-Imada action: 




where the only difference with respect to Spa is in the 
potential energy, which now reads: 

V ef ?(x)=V(x) + ±T(At) 2 (j^j • (18) 

As for the MC move choice, ii), we adopted smarter 
moves, known in the context of classical polymer sim- 
ulations, which reconstruct entire pieces of the polymer, 
instead of a single bead at a time, through the bisection 
algorithm*Z&. This choice guarantees a fast relaxation 
to the instantaneous equilibrium distribution. Moreover, 
we also used global classical moves, in the form of rigid 
motions of the center of mass of the polymer»2J, so as 
to keep a good sampling in the final part of the anneal- 
ing, where the mass m is large (r is small) and quan- 
tum moves are suppressed. In the next Section, we will 
present the annealing results obtained with such an al- 
legedly state-of-the-art PIMC. 

A. PIMC QA Implementation I: IV-order action 
and bisection moves 

Quantum annealing is performed by decreasing the 
inverse-mass parameter, T(t) = S 2 /[2m(t)], appearing in 
Eq. ^| linearly to zero in a time r, T(t) = Tq (1 — t/r) 
(it is understood here that both t and the total annealing 
time r are MC times, i.e., measured in units of MC steps, 
each MC step consisting of one bisection move plus one 
global move). The initial condition is set to To = 0.5 
(as in I). For every value of the annealing time r we cal- 
culated the relevant averages by repeating several times 
the same annealing experiment, starting from different 
randomly distributed initial conditions— 

In Fig. [SI we plot the final PIMC-QA residual energy 
obtained for both potential choices, asymmetric (Top) 
and symmetric (Bottom) (with the same parameters used 
in I and in Sec.[nj, for a fixed temperature T/Vo = 0.03. 
The statistical errors are evaluated by 10 3 repetitions of 
each annealing run for t < 10 5 , and 10 2 repetitions for 
t > 10 5 . We employed P = 160 with an I — 5 bisection 
level (i.e., moving polymer pieces containing 2 5 — 1 = 31 
beads 2 ^, for r < 10 5 , and P = 20 with I — 2 (moving 
2 2 - 1 = 3 beads) for r > 10 5 . We stress that the two 
series of data match perfectly, but obviously the use of 
a smaller P allows us to achieve larger annealing times 
t, which would be otherwise too computationally heavy. 
(An analysis of the convergence behavior of the algorithm 
as a function of both V and P can be found in Ref. 0-) 
In the asymmetric potential case (Top panel of Fig. [SJ 
which, we recall, presents a level crossing and a Landau- 
Zener (LZ) transition (see I) -, the residual energy e res (r) 
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FIG. 8: (Top) PIMC-QA residual energy for the asymmetric 
potential, using a fourth-order action and the bisection algo- 
rithm. The dashed line indicates the thermal equipartition 
limit T/2, for T/Vq = 0.03. As a reference, the results ob- 
tained in Ref. [l^ by exact integration of the imaginary-time 
Schrodinger equation are reported. (Bottom) Same as above, 
for a symmetric potential. 



as a function of the annealing time r gets stuck around 
e res = 0.2. Since this energy is comparable to the energy 
of the metastable state, V(x+), and definitely larger than 
the thermal limit T/2 — 0.015Vb, we see that the al- 
gorithm failed to follow adiabatically the ground state. 
For comparison, we also reported in Fig. |S| the resid- 
ual energy data obtained by the exact imaginary-time 
Schrodinger annealing (reported in I), even if the time- 
scales of the two algorithms are definitely different and 
unrelated. Considering, on the other hand, the symmet- 
ric potential case, bottom panel of Fig. [SJ we notice that 
the algorithm succeeds in reaching the thermal limit T/2 
in a reasonable amount of MC steps. 

From this first comparison between PIMC-QA and ex- 
act Schrodinger QA, we appreciate that the LZ transi- 
tion leads to a severe slowing down of the Monte Carlo 
algorithm, the actual tunneling event being essentially 
missed by the PIMC algorithm. This difficulty is also 
present in static simulations performed in the neighbor- 



hood of the LZ crossing (occurring at T ~ 0.038), where 
we find (results not shown) a dangerous loss of adiabatic- 
ity which calls for an improved sampling of the action: 
We will show how this improved sampling is achieved by 
the introduction of instanton moves. 



B. PIMC Implementation II: Adding the Instanton 

move 

In the previous section we tested an allegedly state- 
of-the-art PIMC-QA algorithm for the very simple prob- 
lem of a particle in a double-well potential, with disap- 
pointing results. The problem is a sampling crisis: Our 
action is accurate, but its sampling misses the tunneling 
events, which is catastrophic when a Landau-Zener cross- 
ing occurs (asymmetric potential case). The well-known 
cure for this kind of sampling problem, for the case of 
a perfectly symmetric double-well, V(x) = Vq(x 2 — a 2 ), 
is the introduction of the so-called instanton move (see 
Ref. l35h . In a nutshell, an instanton is defined as a so- 
lution of the classical equation of motion in the inverted 
potential, which goes from one minimum to the other. 
Moreover, the time-reversed path (anti-instanton) is also 
a solution. Instanton solutions are given by: 



x c i (r) = ± tanh 



Wins (l" — To) 



(19) 



where the ± signs denote instanton and anti-instanton, 



respectively, 



(8TVo)/a is the instanton fre- 



quency, and To a configuration parameter (the instanton 
center). Strictly speaking, the whole classical trajectory 
takes an infinite (real) time, while the barrier overcoming 
is a very fast process (whence the name instanton). In 
our implementation, we made use of the imaginary-time 
version of an instanton/anti-instanton pair as a proposal 
MC mov o 34 ! 35 . In practice, the instanton move proposes 
to a subset of the Trotter's slices an excursion from one 
minimum to the other. The move will be accepted ac- 
cording to the usual Metropolis criterion, with an energy 
competition between the possible gain in potential energy 
compensated by the increase in kinetic energy due to the 
spring stretchings (at the positions of the instanton and 
of the anti-instanton) . Obviously, the instanton move will 
be not effective in the classical limit (small T, final part 
of the annealing), since it describes an inherently quan- 
tum effect: The tunneling process between the two wells. 
In particular, the instanton frequency, LOi ns (defined in 
Eg. I19|) . decreases with T, and when l/u>i ns > (3/2, the 
tunneling time will exceed the thermal cut-off. As a con- 
sequence, the instanton/anti-instanton pair will reduce 
to a useless flat configuration over the interval [0, h(3]. 

Before showing the results, we need to stress an im- 
portant point: The instanton move is available only for 
potentials which are small deformations of a perfectly 
symmetric double-well potential - as our symmetric and 
asymmetric cases are, see Sec. [H] - and is not the gen- 
eral key for solving sampling problems (ergodicity break- 
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ing) for a generic potential, whose landscape is gener- 
ally poorly known. In essence, we are playing here, for 
demonstration purposes, an unfair game, using a detailed 
information on the potential landscape in order to cor- 
rectly implement the tunneling dynamics in our PIMC. 

Fig. [5] shows the PIMC-QA residual energy results 
when the instanton move is introduced, for the case of 
asymmetric (Top) and symmetric potential (Bottom) . As 
in the previous section, we employed P = 160 Trotter's 
slices with I = 5 bisection level, for r < 10 5 , and P = 20 
with I = 2 bisection level for r > 10 5 . Statistical errors 
are evaluated with 10 3 repetitions of every annealing run. 
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the partial acceptance ratio of the instanton moves is 
very low - around 1% -, despite the fact that this pro- 
posal move was "tailored" on the potential under investi- 
gation. Nevertheless, the instanton move quantitatively 
changes the PIMC-QA performance. This is perhaps 
catch a generic feature of systems with barriers, namely 
the crucial importance of rare events - tunneling in this 
case. 

Unfortunately, as said, the instanton "recipe" , natural 
as it is for a double-well potential, cannot be generalized 
to generic landscapes, including more complicated poten- 
tials, not to mention actual combinatorial optimization 
problem. Nevertheless, the results obtained are instruc- 
tive in many respects: They show clearly an important 
limitation of the PIMC approach, while demonstrating, 
once more, that a smart choice of the proposal move - 
even in the case of PIMC-QA - leads to immediate im- 
provement of the annealing performance. 



C. PIMC Implementation III: Choosing the right 
kinetic energy. The Lorentzian move. 

It is by now a sort of Leitmotiv of the paper, that the 
choice of the MC move strongly influences the dynamics, 
and hence the annealing behavior. We recall, in par- 
ticular, see Sec. [HJ that in the classical annealing case 
the winning choice was given by a Lorentzian distributed 
proposal move, which sometimes provides very long dis- 
placements. 

In the present PIMC context, however, the non- 
relativistic form of the kinetic energy part of the Hamil- 
tonian Eq. 1101 forces, in a sense, the choice of Gaussian 
distributed moves, because the free propagator of the 
non-relativistic kinetic energy is a Gaussian, see Eq. ^j] 
The whole bisection algorithm makes strong usaS of the 
Gaussian nature of the free propagator. 

It is natural to ask what would be the QA behavior 
if, instead of the standard non-relativistic kinetic energy 
used so far, we employ for example a relativistic Hamil- 
tonian of the form 



FIG. 9: (Top and Bottom) Same as in Fig. [H] with the in- 
stanton move allowed. As a reference the results obtained 
without instanton move are still reported. The data fit for 
the asymmetric potential case is discussed in the text. 

We note how the introduction of the instanton move 
causes a visible improvement of the residual energy 
slope for the asymmetric potential case: The asymp- 
totic power-law behavior is now quite evident, e res (r) cx 
T -n QA ^ a ^hough with an exponent which is only £Iqa = 
0.19 (smaller than in the Schrodinger integration case, 
where £Iqa = 1/3, see I). For the symmetric potential 
case, the instanton move leads to a faster convergence 
to the thermal limit, but does not really give an over- 
whelming qualitative change. In particular we note that 



H(t) = T(t) \p\ + V(x) 



(20) 



One immediate consequence of this choice is that the free 
Gaussian propagator appearing in Eq.^jjis now replaced 
by a Lorentzian: 



(xi\e 



—A tr\p\ 



Atr 



i 



(21) 



With a certain effort, we have succeeded in generalizing 
the bisection algorithm to implement in a very effective 
way the dynamics provided by this Lorentzian propaga- 
tor (see appendix IbI and Ref. 34 for more details). 

We present here the results of a PIMC-QA approach 
which implements the relativistic kinetic energy, through 
a bisection algorithm with such a "Lorentzian move" . In 
Fig. EH we report the residual energy results for the case 
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of the asymmetric (Top) and symmetric potential (Bot- 
tom). The results are obtained using P — 40 Trotter's 
slices and I = 2 bisection steps. Averages and statisti- 
cal errors are calculated, as usual, with 10 3 repetitions 
of every annealing run. It is clearly seen that, as in the 
classical case, the Lorentzian move - alias, in the present 
context, the relativistic kinetic energy - greatly acceler- 
ates the QA behavior for the difficult asymmetric poten- 
tial case. As usual, the "simpler" symmetric potential 
case does not show a qualitative change. 
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FIG. 10: Same as in Fig. [HI but for a QA based on the rela- 
tivistic kinetic energy of Eq. 1201 implemented via a bisection 
algorithm adapted to Lorentzian moves. As a reference, the 
results obtained by exact integration of the imaginary-time 
Schrodfnger and by the previous Gaussian-based PIMC-QA 
with and without instanton move are still reported. 



transition matrix, thus modifying in a strong way 
the annealing properties. In some sense, the effec- 
tive landscape of the problem is really determined 
not only by the actual potential landscape, but 
also by the choses Markov transition matrix, the 
latter determining the neighborhood of a given 
configuration^. In our double-well example a 
Metropolis MC with Lorentzian-distributed pro- 
posed moves ended up showing the best annealing 
performances, overperforming not only the other 
classical MC choices, but also any Path-Integral 
Monte Carlo (PIMC) algorithm which we were 
able to implement. 

2. Tunneling and sampling difficulties of PIMC. A 

tunneling event, and its associated Landau-Zener 
crossing, can cause severe difficulties to a "state- 
of-the-art" PIMC-QA algorithm. The difficulty 
is associated to a poor sampling of the action 
which misses the rare, but all- important, tunneling 
events. This problem can be cured by the ad-hoc 
introduction of "instanton moves", but this cure, 
while instructive, uses detailed information on the 
potential landscape that is generally not available. 
The generalization of these "instanton moves" to 
more complicated potentials, let alone to generic 
hard optimization problem, seems basically im- 
possible, for it would require, among other things, 
a knowledge of the location of the minima were 
tunneling occurs. 



3. Other limitations of PIMC. PIMC-QA suffers from a 
certain number of other limitations. First of all, it 
works with a finite temperature T, and this sets up 
a thermal energy lower bound below which we can- 
not possibly go (this limit was particularly clear in 
our double- well example). Second, a large number 
of Trotter slices P can cause additional sampling 
problems in that an effective uncorrelation of the 
configurations becomes harder and harder, even if a 
multi-step bisection algorithm is employed. More- 
over, the Trotter break-up itself, see Eq. 1131 can 
cause difficulties whenever, for a generic kinetic en- 
ergy T, the form of the free propagator, general- 
izing the simple Gaussian result of Eq. 1161 is not 
known analytically. This was indeed the difficulty 
met in a PIMC-QA study of the Traveling Sales- 
man Problcmiiii. 



IV. SUMMARY AND DISCUSSION 

Summarizing, we would like to briefly stress some of 
the major conclusions reached in this paper. 

1. Role of move choice in Classical Annealing. The 

choice of the Monte Carlo (MC) moves strongly in- 
fluences the spectral properties of the MC Markov 



4. Role of kinetic energy in QA. The choice of the ki- 
netic energy is clearly all important in QA: 
Sec. IIII CI illustrating the improvements in PIMC- 
QA upon using a relativistic kinetic energy, is par- 
ticularly instructive. In our one-dimensional exam- 
ple, this choice corresponds to a free propagator of 
Lorentzian form, which we implemented by an ap- 
propriately modified bisection scheme (see Adp. IB|| . 
The convenience of such a scheme in more general 
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instances is not a priori obvious and further tests 
on higher-dimensional problems would be needed. 

In view of the previous points, it is fair to say that 
exploring alternative Quantum Monte Carlo schemes for 
performing stochastic implementations of Quantum An- 
nealing remains an important and timely issue for future 
studies in this field. One such scheme, which in principle 
does not suffer from many of the limitations of PIMC, is 
the Green's Function Monte Carlo scheme. Test applica- 
tions of a GFMC-QA strategy to the Ising spin glass are 
currently under way and will be reported shortly^. 
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APPENDIX A: A MODEL FOR THE LOW-LYING 
SPECTRUM OF THE METROPOLIS 
TRANSITION OPERATOR 

In this appendix we shall use analytical tools to explain 
some of the most relevant features of the Markov transi- 
tion matrix spectra presented in Sec. Ill Al The question 
we want to address has to do with the low-lying spectrum 
of a Master equation of the form 

P n+ l(x') = £ W(x'\x) P n (x) , (Al) 

X 

where W(x'|x) is the transition operator associated to a 
certain choice of trial move in a Metropolis MC scheme. 

As we have seen previously (see Sec. Ill A|) . the Box 
proposed move case admits a simple geometrical inter- 
pretation. For the general case, which is much harder 
to treat, we shall sketch here the main steps of a more 



systematic procedure used to extract information on the 
gap of W, referring the reader to Ref. for the many 
technical details. 

As explained in Sec. Ill Al it is possible to find a com- 
plete set of eigenvectors and eigenvalues of the matrix 
W(i\j), defined through the short-time expansion of this 
equation: W = e Atff iS In particular, one can find a lin- 
ear transformation from the (generally) non-symmetric 
operator —IF to a symmetric one, H, which preserves 
its spectral properties^, and is easier to diagonalize in 
a standard way. Moreover, one can apply the same lin- 
ear transformation to the original transition matrix, W; 
this leads to the expansion of the master equation so- 
lution reported in Eq. [8] having taken the appropriate 
time-continuum limit. However, in order to find the W 
eigenvectors and eigenvalues (and in particular its spec- 
tral gap, A, defined in Sec. Ill Al) . it is better to keep At 
finite and consider its simmetrized version, that turns out 
to be the exponential operator e~ AtH , where the appro- 
priate At is uniquely determined by the move range a. 
We notice that the spectrum of —W, {Ai}, introduced in 
Sec. Ill Al is simply related to the spectrum of the original 
transition matrix W, {A^}, by the equation A; = e~ At Ai . 
In the following, we shall refer only to the gap of W. 

The potential V(x) will be always assumed to have the 
usual double-well form. At very small temperatures, the 
equilibrium distribution P eq (x) is concentrated around 
the local minima at x = x± x & . Therefore, there is def- 
initely a regime of parameters where a two-level system 
approximation for W must hold. By exploiting this ap- 
proximation, one can show (see Ref. l34f) that the spectral 
gap A = 1 — Ai of the transition operator W can be ex- 
pressed, at small temperature T, in terms of a well-to- well 
propagator of the evolution operator of H as follows: 

A~ (x^\e- Atn \x + )e' L ^ ± , (A2) 

where V± = V(x±). The crucial quantity we need, there- 
fore, is the well-to- well propagator (x-\e~ AtH \x+). An 
expression for this propagator can be obtained by per- 
forming a Trotter decomposition and writing a Feynman 
path-integral^. The final result can be cast - by an ap- 
propriate saddle-point approximation - in the form: 



(x\e' AtH \x, 



V+);At)e 



-A v— 



P(d(V - V+)) , 



(A3) 



where d(V — V+) is the function introduced in Sec. Ill Al 
(see Fig. 0J), and K\ accounts for the Gaussian fluctu- 
ations around the saddle-point classical solution of the 
Feynman integral, which do not modify the important 
exponentially activated behavior provided by the term 



_ i (y v + +v - \ 
e T V 2 I . P(d(V - V+)), finally, is just the pro- 
posal move distribution: P(d(V — V+)) = P (\xi — X2I) = 
P{xi\x 2 ) (see Fig. El. 

This rather simple equation for the well-to- well propa- 
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gator can be further treated by another standard saddle- 
point approximation of the relevant one-dimensional in- 
tegral. We shall discuss here the cases we are interested 
in: Box, Gaussian, and Lorentzian. 

Box move: P was defined in Eq.Q] We recall that the 
function d(V — V+), introduced in Sec. Ill Al is a monoton- 
ically decreasing function of its argument, which exhibits 
an infinite first order derivative for V = Vo- There is 
no true saddle-point for such a move, since the functions 
involved in the exponential part of Eq. IA3I are all mono- 
tonic. On the other hand, - in the case of the Box move 
- the extremes of integration in Eq. IA3I can be taken as 
[V+, d^ 1 (a) + V + ] instead of [V+, Vo], since its transition 
matrix, P(d(V — V+)), is zero outside the former interval. 
Therefore, the largest contribution to the whole integral 
is due either to V = Vo or V = d -1 (er) + V+, which are 
the two integration extremes. It turns out that the latter 
choice guarantees the largest exponential, and is there- 
fore the right one. As a consequence, we can write down 
that: 



(x\e- AtH \x, 



cx e 



<rV)+- 



(A4) 



This transition amplitude is Arrhenius-like, and so does 
the gap function (see Eq. IA2|) which reads: 



A(T,a) cx e 



(A5) 



where B c s{cr) = d^ 1 (a), as anticipated in Sec. Ill Al from 
geometrical considerations. This is indeed the behavior 
seen in Fig. |3| 

Gaussian move: P was defined in Eq.[5j In this case, 
a true saddle-point is possible if the following equation 
has a solution: 



1 

T 



d(V-V+) d 
a 2 dV 



d(V-V+) = 



(A6) 



Since < and d > 0, then, in the limit T — + we have 
that either V = Vo or V = V+ are solutions of Eq. IA6I 
The second choice guarantees the largest contribution to 
the integral in Eq. IA3I and must be therefore selected. 
The resulting expression for the transition amplitude is: 



other hand, for large values of a we see from Fig. that 
the data deviate from the behavior predicted by Eq. IA8I 
This is the consequence of another level crossing. 

Lorentzian move: P was defined in Eq. To cut 
a long story short (calculations proceed similarly to the 
Box and Gaussian case), we write down the result. The 
gap value is: 



A(T,a) cx a 



(A9) 



and also in this case the model agrees with the data, see 
Fig. [2 within a multiplicative fitting coefficient. 



APPENDIX B: THE BISECTION ALGORITHM 
WITH A LORENTZIAN MOVE 

In this appendix we want to show how to approach, 
within a Path-Integral Monte Carlo, a Hamiltonian with 
relativistic kinetic energy: 



H = T(t) \ P \ + V(x) . 



(Bl) 



The kinetic operator T(t) \p\ is singular in a real space 
representation, but the corresponding density matrix op- 
erator can be written in a simple closed form: 



p K (x',x;At) 



1 



r At 



7T r 2 At 2 + (x' - xY 



(B2) 



This is a Lorentzian (or Cauchy) distribution, while the 
kinetic density matrix operator considered in Sec. II I II 
Eq. 1161 was a Gaussian one. 

It is possible to obtain a generalization of the Levy 
constructional for this kind of Lorentzian move, observ- 
ing that both Gaussian and Lorentzian distribution are 
stable under convolution (a property shared by all the so- 
called Levy distribution*^). The main idea is to obtain a 
closed form of the constrained kinetic operator: 



rp^i), i\ _ PK {x i+ i,x't, At) p K (x' l ,x l - 1 ;At) 

1 K \ X i) ~ 



p K (x i+ i,Xi-i;2At) 



(B3) 



{x. 



\e~ AtH \x, 



(A7) which is the probability of moving the z-th bead to 



and the corresponding gap value is: 



A(T, a) cx e" 



(A8) 



Fig. [3 shows that such a dependence perfectly fits the 
simulation data. We emphasize that even the coeffi- 
cient within the exponential part of the fitting function 
matches the theory (we recall that d(Q) — a cr « 1.75): 
Only an overall prefactor is used as fitting parameter. 

As observed in Sec. Ill Al the previous equation will be 
no longer valid for very small values of the temperature T, 
where an (avoided) crossing between the first eigenvalue 
Ai < 1 and the second eigenvalue A2 will occur. On the 



keeping fixed the (i + l)-th and the (i — l)-th beads in 
Xi + \ and Xi-\, respectively. We recall that this object 
is the building block of the bisection algorithm 27 . In the 
Lorentzian case, it reads: 



Y At 

2tt 



\ TAt 



t x i+1 -x 

\ r At 



/ x[— Xi — 

\ T At 



(B4) 

Unfortunately, this is no longer a Lorentzian distribu- 
tion, and therefore it is not trivial to sample it (for in- 
stance, by making use of the usual transformation tech- 
nique; see Ref . l4Gf) . However, after some algebra, and a 
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smart change of variables (taking due care of the Jaco- 
bian of such a transformation, which changes T onto T), 
we obtain a simpler form: 



4%) = £ 



1 



{y + af + l {y-af + l 



where a 



(x 



i+l 



x i - 1 )/(2TAt), c 



(B5) 

(Xi+l + 



x i - 1 )/(2TAt), and 



TAt 



(B6) 



Finally, with a bit of more algebra, we find that: 



{y + af + l {y-af + l 



{y + af + l (y-af 



(B7) 



w 1 ( v ) 



W 2 ( V ) 



and also that Tjf\y) < 2 W\(y). An efficient strategy to 

sample 2^ (y) is, therefore, to sample W\{y) and then to 
make use of the usual rejection technique (see Ref. Eof) . 
The primitive of Wi{y) can be easily computed. It reads: 

Mi{u) = [ &yWi{y) 

J — oo 

= JLtan-if 2 ^\ + 1 -. (B8) 

Inverting M\{u), and making use of an equidistributed 
random number u G (0, 1], we find that: 




Provided v! S (0, 1] - another equidistributed random 
number - we shall accept the y obtained above according 
to the condition: 

2W 1 (y)u' < W 1 (y)-W 2 (y) . (Bll) 

It results from standard considerations (see Ref. EfTl that 
the average acceptance of this method will be around 
50%. (This means that in order to generate N random 

numbers distributed according to T^\x' i ), one has to try 
on average 2 N times.) Finally, the original middle-point 



value x\ appearing in T K {x[) is obtained by inverting 
the original change of variables in Eq. (|B6(1 . 

x' l =TAt{y + c). (B12) 

This sampling method is a bit cumbersome, but it can be 
easily implemented as a computer algorithm. Its general- 
ization, namely {x[), which represents the probability 
of moving the i-th bead to x\ having fixed the {i + 2 i ~ 1 )- 
th and the (i — 2'~ 1 )-th beads in x i+2 i-i and x i _ 2 i-i, 
can be easily derived (it is enough to change At with 
2 l ~ 1 At). This is what one needs in order to implement 
the bisection scheme^. 

In this way, a PIMC algorithm using such a bisec- 
tion sampling scheme implements correctly a relativis- 
tic choice of Hamiltonian as in Eq. IB1I Moreover, since 
[V, [\p\,V]\ = 0, for such a relativistic-PIMC the sim- 
ple primitive approximation is already exact to fourth- 
order—. We finally report the virial-centroid estimators 
for the kinetic an potential energy for the relativistic- 
PIMC case: 

KpA ~ p + P 2^^- x > dXl 

1 p 

8=1 

where x = 4 Y^i=o x i IS ^he centroid coordinate. The 
way to derive them is completely simil ar to the non- 
relativistic case considered in Refs. I32I41L 
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